% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Replication "Deconstructing the Yield Curve"
% Crump and Gospodinov (2024)
% Date: 26-JUL-2024
% Code to produce Figure 1
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all;
close all;
clc
addpath('functions/');
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Load GSW bond data and produce chart
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
o.dataSet               = 'gsw';
o.freq                  = 'q';
p.N                     = 40;
o.startDate             = 197201;
o.endDate               = 202204;
%
dataAll = loadBondData(p.N);
dataTag = ['dataAll.', o.dataSet, '.', o.freq];
for i = {'prc', 'yld', 'ret', 'xret', 'fwd', 'dret'}, eval([i{1}, 'Raw = ', dataTag, '.', i{1}, ';']);  end
data.mats = (1:p.N);
datesBonds = eval([dataTag, '.dates']);
data.period = eval([dataTag, '.period']);
datesCommon = datesBonds;
%
if ~isempty(o.startDate), datesCommon = datesCommon(datesCommon >= o.startDate); end
if ~isempty(o.endDate), datesCommon = datesCommon(datesCommon <= o.endDate); end
%
dates = datesCommon;
idxBonds = ismember(datesBonds,dates);
%
data.prc = prcRaw(idxBonds,1:p.N);
data.yld = yldRaw(idxBonds,1:p.N);
data.fwd = fwdRaw(idxBonds,1:p.N);
data.dret = dretRaw(idxBonds,1:p.N);
data.ret = retRaw(idxBonds,1:p.N);
data.xret = xretRaw(idxBonds,1:p.N);
data.T = size(data.prc,1);
%
lambda1 = OrderedPCA(data.yld,'corr');
%
figName = figure(1);
plot(lambda1(:,1:3), 'LineWidth', 1.25);
xlim([1 40]);
legend({'PC1','PC2','PC3'},'Interpreter','latex','location','NE', 'Orientation', 'Horizontal')
xlabel('Maturity (Quarters)', 'Interpreter','Latex')
set(gca,'fontsize', 16)
print(figName, '-dpdf', '../output/GSWyld_PCs'); 

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Load oil futures and produce chart
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Monthly data of oil futures for the period Dec. 1999 - Dec. 2016 (1- to
% 24-month maturities)
load ('../input/oilfutures.txt');
year1 = oilfutures(:,1);
month1 = oilfutures(:,2);
loil = log(oilfutures(:,3:end));
lambda2 = OrderedPCA(loil,'corr');
%
figName = figure(2);
plot(lambda2(:,1:3), 'LineWidth', 1.25);
xlim([1 12]);
legend({'PC1','PC2','PC3'},'Interpreter','Latex','location','NE', 'Orientation', 'Horizontal')
xlabel('Maturity (Months)', 'Interpreter','Latex')
set(gca,'fontsize', 16)
print(figName, '-dpdf', '../output/OilFutures_PCs');

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Load S&P 500 options returns and produce chart
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
data_in = '../input/He_Kelly_Manela_Factors_And_Test_Assets.csv';
hkm = importdata(data_in,',');
names = hkm.textdata(1,2:end);
OPTS = hkm.data(66:168,63:80);
OPTS = OPTS(:,10:18);
lambda3 = OrderedPCA(OPTS,'corr');
%
figName = figure(3);
plot(lambda3(:,1:3), 'LineWidth', 1.25);
legend({'PC1','PC2','PC3'},'Interpreter','latex','location','NE', 'Orientation', 'Horizontal')
xlabel('Maturity (Months)', 'Interpreter','Latex')
set(gca,'fontsize', 16)
print(figName, '-dpdf', '../output/SP500Options_PCs');

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Load weather data and produce chart
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
load '../input/weather.txt';
year = (1880:1:2017)';
weather = [weather; zeros(2,2)];
T = size(weather,1);
weath = zeros(length(year),12);
for i=1:size(weath,1)
    weath(i,:) = weather((i-1)*12+1:(i-1)*12+12,2)';
end
lambda4 = pca(zscore(weath));
%
figName = figure(4);
plot([lambda4(:,1:3)], 'LineWidth', 1.25);
xlim([1 12]);
set(gca,'fontsize', 16)
set(gca,'XTick',[1 2 3 4 5 6 7 8 9 10 11 12])
set(gca,'XTickLabel',{'Jan' 'Feb' 'Mar' 'Apr' 'May' 'Jun' 'Jul' 'Aug' 'Sep' 'Oct' 'Nov' 'Dec'},'FontSize',14)
legend({'PC1','PC2','PC3'},'Interpreter','latex', 'location','NE', 'Orientation', 'Horizontal','fontsize', 16)
print(figName, '-dpdf', '../output/Weather_PCs');


